Numerical Treatment of Polar CoordinateSingularities

The treatment of the geometrical singularity in cylindrical and spherical coordinates has for many years been a difficulty in the development of accurate finite difference (FD) and pseudo-spectral (PS) schemes. A variety of numerical procedures for dealing with the singularity have been suggested. For comparative purposes, some of these are discussed in the next sections, but the reader is referred to several books and review papers [3, 7, 10] for more detailed references. Generally, methods discussed in the literature use pole equations, which are akin to boundary conditions to be applied at the singular point. The treatment of the pole as a computational boundary can lead to numerical difficulties. These include the necessity of special boundary closures for FD schemes (e.g., [11]), undesirable clustering of grid points in PS schemes (e.g., [12]), and, in FD schemes, the generation of spurious waves which oscillate from grid point to grid point (so-called two-delta or sawtooth waves, see [4, 26]). In the present paper we investigate a method for treating the coordinate singularity whereby singular coordinates are redefined so that data are differentiated smoothly through the pole, and we avoid placing a grid point directly at the pole. This eliminates the need for any pole equation. Despite the simplicity of the present technique, it appears to be an effective and systematic way to treat many scalar and vector equations in cylindrical and spherical coordinates. A similar technique was used by Merilees [17] for the south and north pole singularities of spherical coordinates but appears not to have been applied more generally. Here we show that the technique leads to excellent results for a number of model problems. The main application we consider is the compressible unsteady Euler and Navier–Stokes equations in cylindrical coordinates (Section 3). For comparison with other methods, we also treat the solution of Bessel’s equation in Section 4. Several